Document last updated 2022-02-03 12:12:27 by Rebecca Shaftel ().

This report summarizes the thermal sensitivity results and draft predictive model for team meeting on February 4, 2022. Goals of the meeting are to discuss model inputs and methods and prediction scenarios.

Objectives:

  1. Build a predictive model for stream thermal sensitivity (TS).
  2. Develop scenarios to understand effects of warming stream temperatures on subsistence salmon populations.

The project mapper is updated with TS by year linked to sites and layers showing a lot of the input data and various hydrologic layers.

1 Thermal sensitivity

DFA output were biased for some years because most of the sites (half or more) came from one watershed with high thermal sensitivity. For some reason, the pattern in these sites all appeared on the trend and their TS was really low.

dfa_dat <- read_csv("DFA/output_8Nov21/GlobalModelEstimates_1trend_DUE.csv")
trend_dat <- read_csv("DFA/output_8Nov21/GlobalTrends_1trend_DUE.csv")

dfa_dat %>% 
  ggplot(aes(TempSens)) +
  geom_freqpoly(aes(color = Region)) +
  facet_wrap(~Year)

trend_dat %>% 
  ggplot(aes(x = DOY, y = Trend)) +
  geom_line() + 
  facet_wrap(~Year)

As an alternative modeling method to estimate stream thermal sensitivity, Tim used the air temperature coefficient in a time series model for each site and summer using an auto-regressive lag of 1 day on the residuals (AR1 model). These results better matched the patterns between air-stream temperatures at the different sites.

mod_dat <- readRDS("data_preparation/final_data/model_data2022-02-02.rds")

mod_dat %>% 
  ggplot(aes(TempSens)) +
  geom_freqpoly(aes(color = Region)) +
  facet_wrap(~Year)

2 Relationships between TS and covariates

panel.cor <- function(x, y, digits = 2, prefix = "", cex.cor, ...)
{
    usr <- par("usr"); on.exit(par(usr))
    par(usr = c(0, 1, 0, 1))
    r <- abs(cor(x, y))
    txt <- format(c(r, 0.123456789), digits = digits)[1]
    txt <- paste0(prefix, txt)
    if(missing(cex.cor)) cex.cor <- 0.8/strwidth(txt)
    text(0.5, 0.5, txt, cex = cex.cor * r)
}

panel.hist <- function(x, ...)
{
    usr <- par("usr"); on.exit(par(usr))
    par(usr = c(usr[1:2], 0, 1.5) )
    h <- hist(x, plot = FALSE)
    breaks <- h$breaks; nB <- length(breaks)
    y <- h$counts; y <- y/max(y)
    rect(breaks[-nB], 0, breaks[-1], y, col = "cyan", ...)
}

mod_dat %>% 
  select(TempSens, log_slope, log_cat_elev, log_cat_slope, 
         wtd_slope_MEAN, snow_ind, wtd_lcld_jd, summer_precip) %>% 
  pairs(upper.panel = panel.cor, diag.panel = panel.hist, lower.panel = panel.smooth)

mod_dat %>% 
  select(TempSens, log_area, dist_coast_km, wtd_north_per,
         asrt_wet, asrt_glac, asrt_lake, glac_10) %>% 
  pairs(upper.panel = panel.cor, diag.panel = panel.hist, lower.panel = panel.smooth)

Exploring interactions between covariates and TS. Starting with different variables and mean watershed slope.

facet_var <- "wtd_slope_MEAN"

mod_dat %>% 
  mutate(fct = cut(get(facet_var), seq(from = min(get(facet_var)) - 0.01, to = max(get(facet_var)), length.out = 5))) %>% 
  ggplot() +
  geom_point(aes(y = TempSens, x = wtd_north_per, color = Region)) +
  geom_smooth(aes(y = TempSens, x = wtd_north_per), method = "lm") +
  facet_wrap(~fct) +
  labs(title = "TS by mean watershed slope and north aspect", subtitle = paste0("Facets are bins of ", facet_var))

facet_var <- "wtd_slope_MEAN"

mod_dat %>% 
  mutate(fct = cut(get(facet_var), seq(from = min(get(facet_var)) - 0.01, to = max(get(facet_var)), length.out = 5))) %>% 
  ggplot() +
  geom_point(aes(y = TempSens, x = asrt_glac, color = Region)) +
  geom_smooth(aes(y = TempSens, x = asrt_glac), method = "lm") +
  facet_wrap(~fct) +
  labs(title = "TS by mean watershed slope and glacier cover", subtitle = paste0("Facets are bins of ", facet_var))

facet_var <- "wtd_slope_MEAN"

mod_dat %>% 
  mutate(fct = cut(get(facet_var), seq(from = min(get(facet_var)) - 0.01, to = max(get(facet_var)), length.out = 5))) %>% 
  ggplot() +
  geom_point(aes(y = TempSens, x = asrt_lake, color = Region)) +
  geom_smooth(aes(y = TempSens, x = asrt_lake), method = "lm") +
  facet_wrap(~fct) +
  labs(title = "TS by mean watershed slope and lake cover", subtitle = paste0("Facets are bins of ", facet_var))

facet_var <- "wtd_slope_MEAN"

mod_dat %>% 
  mutate(fct = cut(get(facet_var), seq(from = min(get(facet_var)) - 0.01, to = max(get(facet_var)), length.out = 5))) %>% 
  ggplot() +
  geom_point(aes(y = TempSens, x = asrt_wet, color = Region)) +
  geom_smooth(aes(y = TempSens, x = asrt_wet), method = "lm") +
  facet_wrap(~fct) +
  labs(title = "TS by mean watershed slope and wetland cover", subtitle = paste0("Facets are bins of ", facet_var))

facet_var <- "wtd_slope_MEAN"

mod_dat %>% 
  mutate(fct = cut(get(facet_var), seq(from = min(get(facet_var)) - 0.01, to = max(get(facet_var)), length.out = 5))) %>% 
  ggplot() +
  geom_point(aes(y = TempSens, x = summer_precip, color = Region)) +
  geom_smooth(aes(y = TempSens, x = summer_precip), method = "lm") +
  facet_wrap(~fct) +
  labs(title = "TS by mean watershed slope and summer precipitation", subtitle = paste0("Facets are bins of ", facet_var))

facet_var <- "wtd_slope_MEAN"

mod_dat %>% 
  mutate(fct = cut(get(facet_var), seq(from = min(get(facet_var)) - 0.01, to = max(get(facet_var)), length.out = 5))) %>% 
  ggplot() +
  geom_point(aes(y = TempSens, x = snow_ind, color = Region)) +
  geom_smooth(aes(y = TempSens, x = snow_ind), method = "lm") +
  facet_wrap(~fct) +
  labs(title = "TS by mean watershed slope and snow index", subtitle = paste0("Facets are bins of ", facet_var))

facet_var <- "wtd_north_per"

mod_dat %>% 
  mutate(fct = cut(get(facet_var), seq(from = min(get(facet_var)) - 0.01, to = max(get(facet_var)), length.out = 5))) %>% 
  ggplot() +
  geom_point(aes(y = TempSens, x = snow_ind, color = Region)) +
  geom_smooth(aes(y = TempSens, x = snow_ind), method = "lm") +
  facet_wrap(~fct) +
  labs(title = "TS by watershed north percent and snow index", subtitle = paste0("Facets are bins of ", facet_var))

facet_var <- "wtd_north_per"

mod_dat %>% 
  mutate(fct = cut(get(facet_var), seq(from = min(get(facet_var)) - 0.01, to = max(get(facet_var)), length.out = 5))) %>% 
  ggplot() +
  geom_point(aes(y = TempSens, x = summer_precip, color = Region)) +
  geom_smooth(aes(y = TempSens, x = summer_precip), method = "lm") +
  facet_wrap(~fct) +
  labs(title = "TS by watershed north percent and summer precipitation", subtitle = paste0("Facets are bins of ", facet_var))

3 Reduce covariates

Reduced the list of covariates based on pairwise correlations and VIF. Also transformed covariates that were heavily right skewed. Arcsine square-root for proportions and log10 for continuous. Glacier cover still heavily right-skewed so added a binary presence-absence for glacier cover as well and checked VIF. Final list of covariates, which are also listed with data sources on google drive sheet:

  • stream order
  • stream gradient
  • distance to coast
  • catchment min, max, and mean slope
  • watershed min, max, and mean slope
  • catchment min, max, and mean elevation
  • watershed min, max, and mean elevation
  • watershed percent north aspect
  • watershed percent wetland cover
  • watershed percent glacier cover
  • watershed percent lake cover
  • watershed area
  • watershed LCLD (last day of longest continuous snow season) by year (2001-2019) - converted to a snow index, the residuals of a model predicting LCLD using mean watershed slope2
  • total summer precipitation at site by year

Removed covariates with pairwise correlations r > 0.7 and VIF > 3. Final list of 13 covariates:

  • log stream gradient
  • log mean catchment elevation
  • log mean catchment slope
  • watershed percent north aspect
  • mean watershed slope
  • log watershed area
  • distance to coast
  • arcsine square-root watershed wetland percent
  • arcsine square-root watershed glacier cover
  • arcsine square-root watershed lake cover
  • snow index
  • summer precipitation
  • glacier presence-absence

4 Select random effects

Created a global model with all 13 covariates plus five interactions from Tim’s Bristol Bay paper and results. Selected best random effects.

Interactions:

  • snow x north
  • slope x north
  • north x precip
  • slope x precip
  • slope x snow

Random effects:

  • none
  • random intercept for Site (multiple TS by year for sites)
  • random intercept for Site within HUC8
  • random intercept for Site within Region (this doesn’t make a lot of sense per mixed effects literature since we only have 5 regions – could include as a fixed effect)

Lowest AIC for HUC8/Site. Predictions can be made at level 1 for HUC8, except for the one or two without data, which we can predict using level 0, or population mean for random intercept.

5 Develop model set

Dredged all subsets of global model to determine variable importance: sum of model weights for all models in which a covariate appears.

lme1_dr <- readRDS("output/lme1m_dredge_results.rds")

#variable importance
lme1_dr %>% 
  as_tibble() %>% 
  mutate(across(asrt_glac:'wtd_north_per:wtd_slope_MEAN', ~ case_when(!is.na(.) ~ weight))) %>% 
  select(asrt_glac:'wtd_north_per:wtd_slope_MEAN') %>% 
  colSums(., na.rm = TRUE) %>% 
  enframe() %>% 
  arrange(desc(value)) %>% 
  ggplot(aes(y = fct_reorder(name, value), x = value)) +
  geom_point() +
  theme_bw() +
  labs(x = "Variable Importance Value", y = "Covariates and Interactions")

Model set of 6 models:

  • global model
  • model with top 12 covariates
  • model with top 8 covariates
  • model with top 4 covariates
  • null model
  • random forest model (I randomly selected only 1 year for each site to avoid pseudo-replicating the spatial covariates)

6 Cross validation

Model selection via cross-validation. Constructing appropriate cross-validation groups to select the optimal model depends on our prediction goals. This paper on cross-validation strategies has some interesting ideas and guidance. I tried three different methods: temporal splits by groups of years, spatial splits by sites in the same HUC10, and leave-one-out.

  • Temporal cross-validation groups were constructed by splitting the years into five groups. For temporal cross-validation, predictions were made for the withheld data using the random intercept for HUC8.
  • Spatial cross-validation was done by constructing 21 groups with 6 HUC10 each that were not adjacent. Predictions were made using the population intercept since these groups overlapped with the HUC8 groupings used for the random intercept.
  • Leave-one-out cross validation withholds one site at a time, builds the model, and predicts for the years of data for that site. I removed any sites that were the only site in a HUC8 in order to predict with the HUC8 random intercepts.

For all cross-validations, the predictions on the testing data were saved and used to calculate RMSE, MAE, and observed versus predicted correlations.

The figure shows the mean RMSE +/- 1SE across cross-validation groups. SE for LOO is missing because groups are only a single site so the RMSE was calculated for all predictions and not by group.

xval_results <- read_rds("output/xval_results.rds") 

xval_results %>% 
  # filter(!xval_type == "LOO") %>% 
  mutate(modelf = factor(model, levels = c("global", "lme_12vars", "lme_8vars", "lme_4vars", "lme_null", "rf"),
                         labels = c("global", "12 vars.", "8 vars.", "4 vars.", "null", "RF"))) %>% 
  ggplot(aes(x = modelf, y = rmse, color = xval_type)) +
  stat_summary(fun = "mean", geom = "point", size = 3) +
  stat_summary(fun.data = "mean_se", geom = "errorbar", width = 0.2) +
  labs(y = "RMSE on Validation Data", x = "Model", color = "Cross-validation Type") +
  theme_bw() +
  theme(legend.position = "bottom", text = element_text(size = 16))

Table summarizing average RMSE, MAE, and correlations for observed versus predicted TS by model.

xval_results %>% 
  group_by(model, xval_type) %>% 
  summarize(mean_rmse = mean(rmse), mean_mae = mean(mae), mean_cor = mean(obs_pred)) %>% 
  kable(digits = 2, col.names = c("Model", "Xval Type", "RMSE", "MAE", "Cor. (r)")) %>% 
  kable_styling()
Model Xval Type RMSE MAE Cor. (r)
global LOO 0.11 0.08 0.48
global spatial 0.11 0.09 0.25
global temporal 0.11 0.08 0.53
lme_12vars LOO 0.11 0.08 0.49
lme_12vars spatial 0.11 0.09 0.29
lme_12vars temporal 0.11 0.08 0.54
lme_4vars LOO 0.11 0.08 0.50
lme_4vars spatial 0.11 0.09 0.27
lme_4vars temporal 0.11 0.08 0.53
lme_8vars LOO 0.11 0.08 0.49
lme_8vars spatial 0.11 0.09 0.28
lme_8vars temporal 0.11 0.08 0.53
lme_null LOO 0.11 0.09 0.45
lme_null spatial 0.12 0.10 NA
lme_null temporal 0.11 0.09 0.48
rf LOO 0.11 0.08 0.52
rf spatial 0.11 0.09 0.27
rf temporal 0.09 0.07 0.71

Plots of observed versus predicted TS by the different models using the LOO cross-validation predictions. The RF model was constructed with just one year for each site as before. Pearson correlations between observed and predicted by region indicate that the models do really poorly for the regions with very little data: Copper, PWS, and Kodiak. We could try region specific models, but we would likely need to consider a smaller set of predictors as these regions don’t have a lot of data.

xval_preds <- readRDS("output/xval_preds.rds")

xval_preds %>% 
  filter(xval_type == "LOO") %>% 
  mutate(modelf = factor(model, levels = c("global", "lme_12vars", "lme_8vars", "lme_4vars", "lme_null", "rf"),
                         labels = c("global", "12 vars.", "8 vars.", "4 vars.", "null", "RF"))) %>% 
  ggplot(aes(x = TempSens, y = preds)) +
  geom_point() +
  geom_abline(aes(intercept = 0, slope = 1)) +
  facet_grid(cols = vars(Region), rows = vars(modelf)) +
  coord_cartesian(xlim = c(-0.2, 1), ylim = c(-0.2,1)) +
  stat_cor(label.x.npc = 0, label.y = c(seq(0.6, 1, 0.1)))

7 CART analysis

Using rpart to create a regression tree with all of the data.

ct_all <- rpart(TempSens ~ log_slope + log_cat_elev + log_cat_slope + wtd_north_per +
                      wtd_slope_MEAN + log_area + dist_coast_km + asrt_wet + 
                      asrt_glac + asrt_lake + snow_ind + summer_precip, 
                    data = mod_dat, method = "anova")

Print the cross validation results to select a complexity parameter for tree-pruning. From plotcp: “A good choice of cp for pruning is often the leftmost value for which the mean lies below the horizontal line.” Select cp for 12 splits.

plotcp(ct_all) 

printcp(ct_all) 
## 
## Regression tree:
## rpart(formula = TempSens ~ log_slope + log_cat_elev + log_cat_slope + 
##     wtd_north_per + wtd_slope_MEAN + log_area + dist_coast_km + 
##     asrt_wet + asrt_glac + asrt_lake + snow_ind + summer_precip, 
##     data = mod_dat, method = "anova")
## 
## Variables actually used in tree construction:
##  [1] asrt_lake      asrt_wet       dist_coast_km  log_area       log_cat_elev  
##  [6] log_cat_slope  snow_ind       summer_precip  wtd_north_per  wtd_slope_MEAN
## 
## Root node error: 26.379/1674 = 0.015758
## 
## n= 1674 
## 
##          CP nsplit rel error  xerror     xstd
## 1  0.120032      0   1.00000 1.00097 0.035900
## 2  0.055597      1   0.87997 0.88256 0.036418
## 3  0.050533      2   0.82437 0.83796 0.036210
## 4  0.037000      3   0.77384 0.79664 0.035599
## 5  0.027847      4   0.73684 0.74375 0.032195
## 6  0.025648      5   0.70899 0.72510 0.032203
## 7  0.020197      7   0.65769 0.68364 0.031625
## 8  0.017486      8   0.63750 0.68468 0.032702
## 9  0.015878      9   0.62001 0.67761 0.032824
## 10 0.011732     10   0.60413 0.66712 0.032316
## 11 0.011346     11   0.59240 0.62882 0.031532
## 12 0.011197     12   0.58106 0.62166 0.031147
## 13 0.010753     19   0.50268 0.62166 0.031147
## 14 0.010000     21   0.48117 0.61138 0.030227
cp_parm <- 0.011197

Prune the tree and print the results.

ct_pr <- prune(ct_all, cp = cp_parm)

# print(ct_pr)

par(xpd = TRUE)
plot(ct_pr, compress = TRUE)
text(ct_pr, use.n = TRUE)

Summary table indicates variable importance and also lists surrogate splits. Interesting that snow index and precipitation are the least important variables.

summary(ct_pr)
## Call:
## rpart(formula = TempSens ~ log_slope + log_cat_elev + log_cat_slope + 
##     wtd_north_per + wtd_slope_MEAN + log_area + dist_coast_km + 
##     asrt_wet + asrt_glac + asrt_lake + snow_ind + summer_precip, 
##     data = mod_dat, method = "anova")
##   n= 1674 
## 
##            CP nsplit rel error    xerror       xstd
## 1  0.12003195      0 1.0000000 1.0009708 0.03590029
## 2  0.05559679      1 0.8799680 0.8825586 0.03641785
## 3  0.05053250      2 0.8243713 0.8379590 0.03620990
## 4  0.03700031      3 0.7738387 0.7966419 0.03559904
## 5  0.02784742      4 0.7368384 0.7437543 0.03219455
## 6  0.02564834      5 0.7089910 0.7250951 0.03220290
## 7  0.02019696      7 0.6576943 0.6836444 0.03162520
## 8  0.01748586      8 0.6374974 0.6846833 0.03270210
## 9  0.01587778      9 0.6200115 0.6776137 0.03282375
## 10 0.01173212     10 0.6041337 0.6671208 0.03231551
## 11 0.01134582     11 0.5924016 0.6288237 0.03153170
## 12 0.01119700     12 0.5810558 0.6216566 0.03114736
## 
## Variable importance
## wtd_slope_MEAN  dist_coast_km       asrt_wet   log_cat_elev       log_area 
##             17             16             14             14              9 
##  log_cat_slope      log_slope      asrt_glac      asrt_lake  wtd_north_per 
##              8              5              5              4              4 
##       snow_ind  summer_precip 
##              3              1 
## 
## Node number 1: 1674 observations,    complexity param=0.120032
##   mean=0.2518006, MSE=0.01575814 
##   left son=2 (885 obs) right son=3 (789 obs)
##   Primary splits:
##       wtd_slope_MEAN < 9.580174    to the right, improve=0.12003200, (0 missing)
##       dist_coast_km  < 83.65199    to the left,  improve=0.08682326, (0 missing)
##       asrt_wet       < 0.279366    to the left,  improve=0.06362352, (0 missing)
##       log_cat_slope  < 1.038689    to the right, improve=0.05420229, (0 missing)
##       log_cat_elev   < 1.912884    to the left,  improve=0.05264191, (0 missing)
##   Surrogate splits:
##       asrt_wet      < 0.1972921   to the left,  agree=0.760, adj=0.490, (0 split)
##       log_cat_slope < 0.7687072   to the right, agree=0.668, adj=0.295, (0 split)
##       asrt_glac     < 0.01209116  to the right, agree=0.662, adj=0.284, (0 split)
##       dist_coast_km < 71.27378    to the left,  agree=0.634, adj=0.223, (0 split)
##       log_slope     < 2.882267    to the right, agree=0.631, adj=0.218, (0 split)
## 
## Node number 2: 885 observations,    complexity param=0.03700031
##   mean=0.210736, MSE=0.01450678 
##   left son=4 (748 obs) right son=5 (137 obs)
##   Primary splits:
##       asrt_wet       < 0.2824055   to the left,  improve=0.07602415, (0 missing)
##       wtd_slope_MEAN < 23.98557    to the right, improve=0.05737510, (0 missing)
##       dist_coast_km  < 123.1085    to the left,  improve=0.04329552, (0 missing)
##       log_area       < 2.096369    to the left,  improve=0.04082977, (0 missing)
##       log_cat_elev   < 2.087551    to the left,  improve=0.03871152, (0 missing)
##   Surrogate splits:
##       dist_coast_km < 205.5014    to the left,  agree=0.855, adj=0.066, (0 split)
##       snow_ind      < -28.03639   to the right, agree=0.855, adj=0.066, (0 split)
##       wtd_north_per < 47.28605    to the left,  agree=0.854, adj=0.058, (0 split)
##       log_area      < 0.5270563   to the right, agree=0.854, adj=0.058, (0 split)
##       log_cat_elev  < 2.842349    to the left,  agree=0.852, adj=0.044, (0 split)
## 
## Node number 3: 789 observations,    complexity param=0.05559679
##   mean=0.2978616, MSE=0.01314866 
##   left son=6 (500 obs) right son=7 (289 obs)
##   Primary splits:
##       log_cat_elev  < 1.912884    to the left,  improve=0.14136820, (0 missing)
##       dist_coast_km < 64.68648    to the left,  improve=0.10994790, (0 missing)
##       log_area      < 0.9796606   to the left,  improve=0.09749469, (0 missing)
##       summer_precip < 1.309348    to the right, improve=0.07405044, (0 missing)
##       asrt_lake     < 0.1675021   to the right, improve=0.06652651, (0 missing)
##   Surrogate splits:
##       dist_coast_km < 77.63445    to the left,  agree=0.897, adj=0.720, (0 split)
##       wtd_north_per < 8.468764    to the right, agree=0.721, adj=0.239, (0 split)
##       asrt_wet      < 0.09142945  to the right, agree=0.702, adj=0.187, (0 split)
##       log_cat_slope < 0.8481061   to the left,  agree=0.677, adj=0.118, (0 split)
##       snow_ind      < 7.872681    to the left,  agree=0.664, adj=0.083, (0 split)
## 
## Node number 4: 748 observations,    complexity param=0.02784742
##   mean=0.1965234, MSE=0.0121109 
##   left son=8 (468 obs) right son=9 (280 obs)
##   Primary splits:
##       log_area       < 2.096369    to the left,  improve=0.08108999, (0 missing)
##       log_cat_elev   < 2.096948    to the left,  improve=0.05740149, (0 missing)
##       wtd_slope_MEAN < 27.78673    to the right, improve=0.05111500, (0 missing)
##       dist_coast_km  < 123.1085    to the left,  improve=0.04888033, (0 missing)
##       wtd_north_per  < 7.072823    to the left,  improve=0.04036520, (0 missing)
##   Surrogate splits:
##       asrt_glac     < 0.01156298  to the left,  agree=0.699, adj=0.196, (0 split)
##       dist_coast_km < 90.30846    to the left,  agree=0.686, adj=0.161, (0 split)
##       log_slope     < 2.289019    to the right, agree=0.678, adj=0.139, (0 split)
##       snow_ind      < 12.74804    to the left,  agree=0.678, adj=0.139, (0 split)
##       asrt_wet      < 0.1076148   to the left,  agree=0.652, adj=0.071, (0 split)
## 
## Node number 5: 137 observations,    complexity param=0.02019696
##   mean=0.2883342, MSE=0.02046354 
##   left son=10 (59 obs) right son=11 (78 obs)
##   Primary splits:
##       wtd_slope_MEAN < 13.90976    to the left,  improve=0.1900400, (0 missing)
##       snow_ind       < -2.388058   to the right, improve=0.1894510, (0 missing)
##       log_cat_elev   < 1.149802    to the right, improve=0.1644225, (0 missing)
##       dist_coast_km  < 20.35841    to the right, improve=0.1464481, (0 missing)
##       wtd_north_per  < 10.67759    to the left,  improve=0.1129494, (0 missing)
##   Surrogate splits:
##       asrt_wet      < 0.3845997   to the right, agree=0.839, adj=0.627, (0 split)
##       log_area      < 0.6163236   to the left,  agree=0.774, adj=0.475, (0 split)
##       wtd_north_per < 16.42505    to the left,  agree=0.766, adj=0.458, (0 split)
##       log_cat_elev  < 1.500962    to the left,  agree=0.723, adj=0.356, (0 split)
##       dist_coast_km < 0.1619185   to the left,  agree=0.672, adj=0.237, (0 split)
## 
## Node number 6: 500 observations,    complexity param=0.02564834
##   mean=0.2650838, MSE=0.01112529 
##   left son=12 (137 obs) right son=13 (363 obs)
##   Primary splits:
##       log_area      < 1.723963    to the left,  improve=0.10713680, (0 missing)
##       summer_precip < 1.061576    to the right, improve=0.09602527, (0 missing)
##       log_cat_elev  < 1.382246    to the left,  improve=0.09384839, (0 missing)
##       asrt_lake     < 0.05347914  to the right, improve=0.07654909, (0 missing)
##       wtd_north_per < 24.39462    to the left,  improve=0.07173712, (0 missing)
##   Surrogate splits:
##       summer_precip < 5.381685    to the right, agree=0.798, adj=0.263, (0 split)
##       log_cat_elev  < 1.328695    to the left,  agree=0.792, adj=0.241, (0 split)
##       log_slope     < 2.883745    to the right, agree=0.786, adj=0.219, (0 split)
##       dist_coast_km < 0.571042    to the left,  agree=0.784, adj=0.212, (0 split)
##       snow_ind      < 14.10572    to the right, agree=0.772, adj=0.168, (0 split)
## 
## Node number 7: 289 observations,    complexity param=0.0505325
##   mean=0.3545707, MSE=0.01157457 
##   left son=14 (33 obs) right son=15 (256 obs)
##   Primary splits:
##       log_cat_elev   < 2.560219    to the right, improve=0.3985000, (0 missing)
##       dist_coast_km  < 133.6306    to the right, improve=0.3935184, (0 missing)
##       wtd_slope_MEAN < 1.027043    to the left,  improve=0.1806864, (0 missing)
##       log_cat_slope  < -0.03784077 to the left,  improve=0.1146985, (0 missing)
##       asrt_lake      < 0.001917876 to the left,  improve=0.1065575, (0 missing)
##   Surrogate splits:
##       dist_coast_km  < 106.8269    to the right, agree=0.983, adj=0.848, (0 split)
##       wtd_slope_MEAN < 1.027043    to the left,  agree=0.920, adj=0.303, (0 split)
##       log_cat_slope  < 1.031943    to the right, agree=0.903, adj=0.152, (0 split)
##       asrt_wet       < 0.007257437 to the left,  agree=0.903, adj=0.152, (0 split)
##       log_area       < 3.541905    to the right, agree=0.896, adj=0.091, (0 split)
## 
## Node number 8: 468 observations
##   mean=0.1722837, MSE=0.01081693 
## 
## Node number 9: 280 observations,    complexity param=0.01587778
##   mean=0.2370384, MSE=0.01165014 
##   left son=18 (66 obs) right son=19 (214 obs)
##   Primary splits:
##       log_cat_slope < 0.4658453   to the left,  improve=0.1283988, (0 missing)
##       asrt_lake     < 0.145426    to the right, improve=0.1233636, (0 missing)
##       log_cat_elev  < 2.10936     to the left,  improve=0.1186668, (0 missing)
##       asrt_glac     < 0.2219403   to the right, improve=0.1172669, (0 missing)
##       log_slope     < 2.601229    to the left,  improve=0.1083309, (0 missing)
##   Surrogate splits:
##       log_area     < 3.882179    to the right, agree=0.839, adj=0.318, (0 split)
##       log_cat_elev < 1.486185    to the left,  agree=0.789, adj=0.106, (0 split)
##       asrt_glac    < 0.5472557   to the right, agree=0.789, adj=0.106, (0 split)
##       asrt_lake    < 0.01812404  to the left,  agree=0.779, adj=0.061, (0 split)
##       asrt_wet     < 0.1379552   to the right, agree=0.775, adj=0.045, (0 split)
## 
## Node number 10: 59 observations,    complexity param=0.01748586
##   mean=0.2166317, MSE=0.01555308 
##   left son=20 (43 obs) right son=21 (16 obs)
##   Primary splits:
##       asrt_lake      < 0.25946     to the left,  improve=0.5026653, (0 missing)
##       wtd_north_per  < 14.58423    to the left,  improve=0.3189109, (0 missing)
##       summer_precip  < 4.889728    to the left,  improve=0.3016230, (0 missing)
##       dist_coast_km  < 13.10007    to the right, improve=0.2904042, (0 missing)
##       wtd_slope_MEAN < 10.74938    to the right, improve=0.2900983, (0 missing)
##   Surrogate splits:
##       wtd_slope_MEAN < 10.74938    to the right, agree=0.898, adj=0.625, (0 split)
##       log_slope      < 0.9154418   to the right, agree=0.864, adj=0.500, (0 split)
##       log_area       < 0.3408574   to the right, agree=0.864, adj=0.500, (0 split)
##       asrt_wet       < 0.3626343   to the right, agree=0.847, adj=0.438, (0 split)
##       wtd_north_per  < 14.58423    to the left,  agree=0.797, adj=0.250, (0 split)
## 
## Node number 11: 78 observations,    complexity param=0.01173212
##   mean=0.3425707, MSE=0.01734737 
##   left son=22 (15 obs) right son=23 (63 obs)
##   Primary splits:
##       snow_ind      < 10.2671     to the right, improve=0.2287222, (0 missing)
##       asrt_wet      < 0.376018    to the left,  improve=0.2109144, (0 missing)
##       log_area      < 0.8209101   to the right, improve=0.1798855, (0 missing)
##       log_cat_elev  < 1.18285     to the right, improve=0.1798855, (0 missing)
##       dist_coast_km < 5.586949    to the right, improve=0.1627627, (0 missing)
##   Surrogate splits:
##       wtd_north_per < 35.90707    to the right, agree=0.859, adj=0.267, (0 split)
##       log_area      < 3.490255    to the right, agree=0.833, adj=0.133, (0 split)
##       asrt_glac     < 0.2351146   to the right, agree=0.833, adj=0.133, (0 split)
## 
## Node number 12: 137 observations
##   mean=0.2088862, MSE=0.008942191 
## 
## Node number 13: 363 observations,    complexity param=0.02564834
##   mean=0.2862934, MSE=0.01030744 
##   left son=26 (276 obs) right son=27 (87 obs)
##   Primary splits:
##       dist_coast_km < 1.44316     to the right, improve=0.20237260, (0 missing)
##       log_cat_slope < 0.6440247   to the left,  improve=0.15974730, (0 missing)
##       asrt_lake     < 0.05670925  to the right, improve=0.13122180, (0 missing)
##       summer_precip < 1.061576    to the right, improve=0.09970441, (0 missing)
##       wtd_north_per < 24.49116    to the left,  improve=0.07031372, (0 missing)
##   Surrogate splits:
##       asrt_lake     < 0.04391394  to the right, agree=0.917, adj=0.655, (0 split)
##       log_cat_slope < 0.6440247   to the left,  agree=0.895, adj=0.563, (0 split)
##       log_slope     < 2.763205    to the left,  agree=0.813, adj=0.218, (0 split)
##       summer_precip < 0.7691304   to the right, agree=0.796, adj=0.149, (0 split)
##       wtd_north_per < 30.37812    to the left,  agree=0.788, adj=0.115, (0 split)
## 
## Node number 14: 33 observations
##   mean=0.1654106, MSE=0.005788662 
## 
## Node number 15: 256 observations
##   mean=0.3789547, MSE=0.007113369 
## 
## Node number 18: 66 observations
##   mean=0.1673949, MSE=0.01090623 
## 
## Node number 19: 214 observations
##   mean=0.2585173, MSE=0.009922364 
## 
## Node number 20: 43 observations
##   mean=0.1626964, MSE=0.004643366 
## 
## Node number 21: 16 observations
##   mean=0.3615829, MSE=0.01604409 
## 
## Node number 22: 15 observations
##   mean=0.2134797, MSE=0.00267514 
## 
## Node number 23: 63 observations
##   mean=0.3733066, MSE=0.01592834 
## 
## Node number 26: 276 observations,    complexity param=0.01134582
##   mean=0.2606511, MSE=0.008576905 
##   left son=52 (41 obs) right son=53 (235 obs)
##   Primary splits:
##       log_cat_elev  < 1.386769    to the left,  improve=0.12643180, (0 missing)
##       asrt_lake     < 0.1657009   to the right, improve=0.12421390, (0 missing)
##       summer_precip < 1.088696    to the right, improve=0.08968576, (0 missing)
##       dist_coast_km < 41.47178    to the left,  improve=0.07918064, (0 missing)
##       log_area      < 3.352934    to the right, improve=0.05877336, (0 missing)
##   Surrogate splits:
##       asrt_lake < 0.3193176   to the right, agree=0.909, adj=0.390, (0 split)
##       log_area  < 3.485453    to the right, agree=0.862, adj=0.073, (0 split)
##       asrt_glac < 0.0481324   to the right, agree=0.862, adj=0.073, (0 split)
## 
## Node number 27: 87 observations
##   mean=0.3676412, MSE=0.007093985 
## 
## Node number 52: 41 observations
##   mean=0.1818131, MSE=0.009323335 
## 
## Node number 53: 235 observations
##   mean=0.2744058, MSE=0.007173091

8 Post-hoc exploration of regional TS

For each region, the number of sites is listed along with a table showing the years for each site, and pairplots of TS against each of the covariates (separated into two plots since there are so many).

Copper River

region_in <- "Copper_River"

print(paste0(region_in, ": ", nrow(mod_dat %>% filter(Region == region_in) %>% distinct(Site)), " sites"))
## [1] "Copper_River: 28 sites"
mod_dat %>% 
  filter(Region == region_in) %>% 
  group_by(Site) %>% 
  summarize(years = paste(Year, collapse = ", ")) %>% 
  kable(longtable = TRUE) %>%
  kable_styling()
Site years
fws_Long Lake Creek 2009, 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018
fws_Tanada Creek 2009, 2010, 2013, 2014, 2015, 2016, 2017, 2019
NPS_Caribou Creek 2008, 2009, 2010, 2011, 2012
NPS_Crystal Creek 2011, 2012, 2013, 2014
NPS_Gilahina River 2008, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019
NPS_Lakina River 2008
NPS_Long Lake Creek 2008, 2011, 2012, 2013, 2014
NPS_Rock Creek WRST 2010, 2011, 2012, 2018
NPS_Rufus Creek 2009, 2010, 2011, 2012, 2014, 2015, 2016, 2017, 2018, 2019
USFS_18 Mile 2019
USFS_18 Mile Middle Fork 2013, 2014, 2017, 2018
USFS_18 Mile West Fork 2013, 2017
USFS_24.9 Mile Creek 2013, 2014, 2016, 2017, 2018
USFS_25 Mile 2011, 2012, 2014, 2015, 2016, 2017, 2018, 2019
USFS_Blackhole Creek 2014, 2015, 2016, 2017, 2018, 2019
USFS_Cabin Lake Outlet 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019
USFS_East Fork 18 Mile 2010, 2011, 2012, 2014, 2015, 2016, 2017, 2018
USFS_Eyak Lake Tributary 2019
USFS_Hook Point 2014, 2015, 2016, 2017, 2018, 2019
USFS_Ibeck Creek-Low 2016, 2017, 2018, 2019
USFS_Ibeck Creek-Lower Side Channel 2013, 2014
USFS_Little Martin River 2011, 2012, 2013, 2014, 2016, 2017, 2018, 2019
USFS_Martin Lake- Inlet 2011, 2012, 2013, 2014, 2016, 2017, 2018, 2019
USFS_Middle Arm Eyak 2012, 2014, 2015, 2016, 2017, 2018
USFS_Power Creek 2010, 2011, 2012, 2013, 2015, 2016, 2017, 2018, 2019
USFS_Salmon Creek 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019
usgs_15200280 2016, 2018, 2019
usgs_15215900 2013, 2014, 2015, 2016, 2017, 2018, 2019
mod_dat %>% 
  filter(Region == region_in) %>% 
  select(TempSens, log_slope, log_cat_elev, log_cat_slope, 
         wtd_slope_MEAN, snow_ind, wtd_lcld_jd, summer_precip) %>% 
  pairs(upper.panel = panel.cor, diag.panel = panel.hist, lower.panel = panel.smooth)

mod_dat %>% 
  filter(Region == region_in) %>% 
  select(TempSens, log_area, dist_coast_km, wtd_north_per,
         asrt_wet, asrt_glac, asrt_lake, glac_10) %>% 
  pairs(upper.panel = panel.cor, diag.panel = panel.hist, lower.panel = panel.smooth)

Kodiak

region_in <- "Kodiak"

print(paste0(region_in, ": ", nrow(mod_dat %>% filter(Region == region_in) %>% distinct(Site)), " sites"))
## [1] "Kodiak: 31 sites"
mod_dat %>% 
  filter(Region == region_in) %>% 
  group_by(Site) %>% 
  summarize(years = paste(Year, collapse = ", ")) %>% 
  kable(longtable = TRUE) %>%
  kable_styling()
Site years
571005154134600 2005, 2006
571221154040300 2005
571343154244900 2006
572656154082400 2005, 2006
fws_kdk_aforv01 2009, 2010, 2011, 2012, 2013, 2014, 2016, 2017, 2018
fws_kdk_busrv01 2009, 2012, 2014, 2015, 2016, 2017
kdk_akacr01 2016, 2017
kdk_akacr02 2018
kdk_ayarv01 2016, 2017, 2018
kdk_ayarv03 2015, 2016, 2017, 2018
kdk_bigcr01 2016, 2018
kdk_busrv01 2016, 2017, 2018
kdk_cancr01 2015, 2016, 2017
kdk_cancr02 2018
kdk_cascr01 2015, 2016, 2017, 2018, 2019
kdk_concr01 2015, 2016, 2017, 2018, 2019
kdk_doscr01 2015, 2016, 2017, 2018, 2019
kdk_doscr02 2015, 2017, 2018, 2019
kdk_eftrv01 2015, 2016, 2017, 2018
kdk_falcr01 2015, 2017
kdk_karrv01 2015, 2016, 2017, 2018, 2019
kdk_karrv02 2016, 2017, 2018
kdk_meacr01 2015, 2016
kdk_olgcr01a 2016, 2017, 2019
kdk_omarv01 2015, 2016
kdk_pincr01 2015, 2016, 2017
kdk_pincr02 2018
kdk_relrv01 2015, 2016, 2017, 2018
kdk_soucr01 2015, 2017, 2018
usgs_15295700 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019
usgs_15297475 2018, 2019
mod_dat %>% 
  filter(Region == region_in) %>% 
  select(TempSens, log_slope, log_cat_elev, log_cat_slope, 
         wtd_slope_MEAN, snow_ind, wtd_lcld_jd, summer_precip) %>% 
  pairs(upper.panel = panel.cor, diag.panel = panel.hist, lower.panel = panel.smooth)

mod_dat %>% 
  filter(Region == region_in) %>% 
  select(TempSens, log_area, dist_coast_km, wtd_north_per,
         asrt_wet, asrt_glac, asrt_lake, glac_10) %>% 
  pairs(upper.panel = panel.cor, diag.panel = panel.hist, lower.panel = panel.smooth)

Prince William Sound

region_in <- "Prince_William_Sound"

print(paste0(region_in, ": ", nrow(mod_dat %>% filter(Region == region_in) %>% distinct(Site)), " sites"))
## [1] "Prince_William_Sound: 20 sites"
mod_dat %>% 
  filter(Region == region_in) %>% 
  group_by(Site) %>% 
  summarize(years = paste(Year, collapse = ", ")) %>% 
  kable(longtable = TRUE) %>%
  kable_styling()
Site years
nv_Eyak_Hartney Creek 2016, 2017
nv_Eyak_Heney Creek 2017, 2018
PWSCC_Erb 2016, 2017, 2018
PWSCC_Gilmour 2016, 2019
PWSCC_Hogan 2018, 2019
PWSCC_Stockdale 2016, 2017, 2018, 2019
USFS_Eagle Creek 2014, 2015, 2016, 2017, 2018, 2019
USFS_ERB Creek 2019
USFS_Hell’s Hole Trib 2013, 2014, 2015, 2016, 2017, 2018, 2019
USFS_Jackpot River 2014, 2015, 2016, 2017
USFS_Koppen Creek 2013, 2014, 2015, 2016, 2017, 2018, 2019
USFS_Olsen Creek 2013, 2014, 2015, 2016, 2017, 2018
USFS_Pigot Bay Spawn Channel 2014, 2015, 2016, 2017
USFS_Rude River SC 2014, 2015, 2017, 2018
USFS_Sheep River 2013, 2014, 2015, 2016, 2017, 2018, 2019
USFS_Shelter Bay Trib 2013, 2014, 2015, 2016, 2017, 2018
USFS_Solf Lake Fish Pass 2012, 2013, 2014, 2015, 2016, 2017, 2018
USFS_Solf Lake Outlet Creek 2019
USFS_Stump Lake Outlet 2014, 2015, 2016, 2017, 2018
usgs_15236900 2016, 2017, 2018, 2019
mod_dat %>% 
  filter(Region == region_in) %>% 
  select(TempSens, log_slope, log_cat_elev, log_cat_slope, 
         wtd_slope_MEAN, snow_ind, wtd_lcld_jd, summer_precip) %>% 
  pairs(upper.panel = panel.cor, diag.panel = panel.hist, lower.panel = panel.smooth)

mod_dat %>% 
  filter(Region == region_in) %>% 
  select(TempSens, log_area, dist_coast_km, wtd_north_per,
         asrt_wet, asrt_glac, asrt_lake, glac_10) %>% 
  pairs(upper.panel = panel.cor, diag.panel = panel.hist, lower.panel = panel.smooth)

Cook Inlet

region_in <- "Cook_Inlet"

print(paste0(region_in, ": ", nrow(mod_dat %>% filter(Region == region_in) %>% distinct(Site)), " sites"))
## [1] "Cook_Inlet: 211 sites"
mod_dat %>% 
  filter(Region == region_in) %>% 
  group_by(Site) %>% 
  summarize(years = paste(Year, collapse = ", ")) %>% 
  kable(longtable = TRUE) %>%
  kable_styling()
Site years
APU10 2015
APU11 2015
APU12 2015
APU5 2015
AR2 2006, 2007
BM1 2015, 2016
BM10 2015, 2016
BM11 2016
BM12 2016
BM13 2015
BM14 2016
BM2 2015, 2016
BM3 2015
BM6 2015
BM7 2016
BM8 2015, 2016
BM9 2015, 2016
CIK_0 2009, 2010, 2011, 2012
CIK_1 2008, 2009, 2010, 2011
CIK_10 2008, 2009, 2010, 2011, 2012, 2019
CIK_11 2008, 2009, 2010, 2011, 2012, 2014, 2015, 2019
CIK_12 2008, 2009, 2010, 2011, 2012, 2014, 2015
CIK_13 2008, 2009, 2010, 2011, 2014, 2015
CIK_14 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019
CIK_15 2004, 2005, 2006, 2007, 2008, 2009, 2011
CIK_16 2009, 2010, 2011, 2012
CIK_17 2008, 2009, 2010, 2011, 2012
CIK_18 2009, 2010, 2011, 2012
CIK_19 2009, 2010, 2011, 2012
CIK_2 2008, 2009, 2010
CIK_20 2008, 2009, 2010, 2011, 2012
CIK_21 2008, 2009, 2010, 2011, 2012
CIK_22 2009, 2010
CIK_23 2008, 2009, 2010, 2011, 2012, 2014, 2015, 2019
CIK_24 2009, 2010, 2011, 2012
CIK_25 2008, 2010, 2012
CIK_26 2009, 2010, 2011, 2012
CIK_27 2008, 2009, 2010, 2011
CIK_28 2009, 2010, 2014, 2016, 2017
CIK_29 2008, 2009, 2010, 2012
CIK_3 2008, 2009, 2010, 2011, 2013, 2014, 2015, 2016, 2017, 2018, 2019
CIK_30 2008, 2009, 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019
CIK_31 2008, 2009, 2010, 2011, 2012
CIK_32 2009, 2010, 2011
CIK_33 2008, 2009, 2010, 2011
CIK_34 2009, 2011
CIK_35 2008, 2009, 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019
CIK_36 2008, 2009, 2010, 2011
CIK_37 2008, 2009, 2010, 2011
CIK_38 2008, 2009, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019
CIK_39 2008, 2009, 2010, 2011, 2012
CIK_4 2008, 2009, 2010, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019
CIK_40 2008, 2009, 2010, 2011, 2012
CIK_41 2008, 2009, 2010, 2012
CIK_42 2009, 2010
CIK_43 2008, 2009, 2011, 2012
CIK_44 2009
CIK_45 2011
CIK_46 2010, 2011
CIK_47 2009, 2010
CIK_48 2011
CIK_5 2008, 2009, 2010, 2011, 2012, 2014, 2015
CIK_6 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2013, 2014, 2015, 2016, 2017, 2018, 2019
CIK_7 2008, 2009, 2010, 2011
CIK_8 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019
CIK_9 2008, 2009, 2010, 2011, 2012, 2014, 2015, 2019
CIK1 2005, 2006, 2007
CIK2 2005
CIK3 2005, 2006, 2007
CIK4 2005, 2006, 2007
CIK6 2005
Deshka 1 Downstream 2017, 2018, 2019
Deshka 1 Tributary 2017, 2018, 2019
Deshka 1 Upstream 2017, 2018
Deshka 2 Downstream 2018, 2019
Deshka 2 Tributary 2017, 2018, 2019
Deshka 2 Upstream 2017, 2018
Deshka 3 Downstream 2017, 2018, 2019
Deshka 3 Tributary 2017, 2018, 2019
Deshka 3 Upstream 2018
Deshka 4 Downstream 2017, 2019
Deshka 4 Tributary 2017, 2018, 2019
Deshka 4 Upstream 2017, 2018, 2019
Deshka 5 Downstream 2018, 2019
Deshka 5 Tributary 2017, 2018, 2019
Deshka 6 Downstream 2017, 2018, 2019
Deshka 6 Tributary 2017, 2018, 2019
Deshka 6 Upstream 2017, 2018, 2019
Kroto 1 Downstream 2017, 2018, 2019
Kroto 1 Tributary 2017, 2018, 2019
Kroto 1 Upstream 2017, 2018, 2019
Kroto 2 Downstream 2017, 2018, 2019
Kroto 2 Tributary 2017, 2018, 2019
Kroto 2 Upstream 2017, 2018, 2019
Kroto 3 Downstream 2017, 2018, 2019
Kroto 3 Tributary 2017, 2018, 2019
Kroto 3 Upstream 2017, 2018, 2019
Kroto 4 Downstream 2017, 2018, 2019
Kroto 4 Tributary 2017, 2018, 2019
Kroto 4 Upstream 2017, 2018, 2019
Kroto 5 Downstream 2017, 2018, 2019
Kroto 6 Downstream 2017, 2018, 2019
Kroto 6 Tributary 2018
Kroto 6 Upstream 2017, 2018, 2019
Kroto 7 Downstream 2017, 2018, 2019
Kroto 7 Tributary 2018, 2019
Kroto 7 Upstream 2017, 2018, 2019
Kroto 8 Downstream 2017, 2018
Kroto 8 Tributary 2017, 2018
Kroto 8 Upstream 2017, 2018, 2019
KWF1 2015, 2016
KWF2 2015, 2016
KWF3 2016
lsarc 2016
lscoh 2016
lslak 2016
lslak7 2016
lslil10 2016
lspap 2016
lsr103 2016
lsr112.2 2016
lsr112.4 2016
lsr28 2016
lsr33.1 2016
lsr33.8 2016
lsr48.1 2016
lsr48.3 2016
lsr63.1 2016
lsr63.4 2016
lsr73 2016
lsr73.1 2016
lsr87 2016
lstrb33 2016
Moose 1 Top of Study 2017, 2018, 2019
Moose 2 Downstream 2017, 2018, 2019
Moose 2 Tributary 2017, 2018, 2019
Moose 2 Upstream 2017, 2018, 2019
Moose 3 Downstream 2017, 2018, 2019
Moose 3 Tributary 2017, 2018
Moose 3 Upstream 2017, 2019
Moose 4 Downstream 2017, 2018, 2019
Moose 4 Tributary 2018, 2019
Moose 4 Upstream 2017, 2018, 2019
Moose 5 Downstream 2017, 2018, 2019
Moose 5 Tributary 2018, 2019
Moose 5 Upstream 2017, 2018
Moose 6 Downstream 2017, 2018, 2019
Moose 6 Tributary 2017, 2018, 2019
Moose 6 Upstream 2017, 2018, 2019
Moose 7 Downstream 2017, 2018, 2019
Moose 7 Tributary 2019
Moose 7 Upstream 2017, 2018, 2019
Moose 8 Downstream 2018, 2019
Moose 8 Tributary 2017, 2018, 2019
Moose 8 Upstream 2017, 2018, 2019
OMCT1 2019
OMCT2 2019
OMCT3 2019
OMCT4 2019
OWL1 2019
PGC1 2019
PGC2 2019
PGCT1 2019
PJ Mid 2012
PMC1 2019
PMC2 2019
PMCT1 2019
PNC1 2019
PNCT1 2019
PNCT2 2019
PNCT4 2019
PWMC1 2019
STAR 171 Lower 2012
STAR 171 Upper 2012
USFS_Bench Creek 2014, 2015, 2016, 2018, 2019
USFS_Center Creek 2015, 2016, 2017, 2018, 2019
USFS_Chickaloon Headwaters 2014, 2016, 2017
USFS_Crescent Creek 2014, 2015, 2016, 2017, 2018
USFS_Daves Creek 2014, 2015, 2016, 2017, 2018, 2019
USFS_Juneau Creek 2014, 2015, 2016, 2017, 2018
USFS_Quartz Creek 2014, 2015, 2016, 2017, 2018, 2019
USFS_Resurrection Creek 2014, 2015, 2016, 2017, 2018
usgs_15238984 2013
usgs_15238986 2011, 2018, 2019
usgs_15239070 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019
usgs_15241600 2001, 2002
usgs_15261000 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019
usgs_15266110 2001
usgs_15266300 2001, 2014, 2015, 2016, 2017, 2018, 2019
usgs_15272380 2003, 2004, 2006, 2007, 2008, 2009, 2010, 2011
usgs_15274000 2001
usgs_15274600 2018, 2019
usgs_15276000 2015, 2016, 2017, 2018, 2019
usgs_15280999 2019
usgs_15283700 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019
usgs_15284000 2015, 2016, 2017, 2018, 2019
usgs_15285000 2010, 2011, 2012
usgs_15290000 2015, 2016, 2017, 2018, 2019
usgs_15291000 2012, 2013, 2014
usgs_15291700 2016
usgs_15292100 2012, 2014
usgs_15292400 2012, 2013, 2014, 2016
usgs_15292700 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019
usgs_15292780 2012, 2013, 2014, 2015
usgs_15292800 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018
usgs_15293200 2013, 2014, 2015, 2016, 2017, 2018
usgs_15293700 2013, 2014, 2015, 2016, 2017, 2018
usgs_15294005 2015, 2016, 2017, 2018, 2019
usgs_15294080 2017, 2018, 2019
usgs_15294345 2013, 2014
usgs_15294350 2013, 2014
mod_dat %>% 
  filter(Region == region_in) %>% 
  select(TempSens, log_slope, log_cat_elev, log_cat_slope, 
         wtd_slope_MEAN, snow_ind, wtd_lcld_jd, summer_precip) %>% 
  pairs(upper.panel = panel.cor, diag.panel = panel.hist, lower.panel = panel.smooth)

mod_dat %>% 
  filter(Region == region_in) %>% 
  select(TempSens, log_area, dist_coast_km, wtd_north_per,
         asrt_wet, asrt_glac, asrt_lake, glac_10) %>% 
  pairs(upper.panel = panel.cor, diag.panel = panel.hist, lower.panel = panel.smooth)

Bristol Bay

region_in <- "Bristol_Bay"

print(paste0(region_in, ": ", nrow(mod_dat %>% filter(Region == region_in) %>% distinct(Site)), " sites"))
## [1] "Bristol_Bay: 118 sites"
mod_dat %>% 
  filter(Region == region_in) %>% 
  group_by(Site) %>% 
  summarize(years = paste(Year, collapse = ", ")) %>% 
  kable(longtable = TRUE) %>%
  kable_styling()
Site years
accs_AKBB-011 2015, 2016, 2017, 2018, 2019
accs_AKBB-020 2016, 2017, 2018, 2019
accs_AKBB-028 2015, 2016, 2017, 2018, 2019
accs_AKBB-029 2015, 2016, 2017, 2018, 2019
accs_AKBB-040 2016, 2017, 2018, 2019
accs_ilbig01 2018
accs_iltaz02 2016, 2017, 2018, 2019
accs_iltnr19 2013, 2014, 2015, 2016, 2017, 2018, 2019
accs_mubon10 2015, 2016, 2019
accs_muekm23 2013, 2014, 2015, 2016, 2017, 2018, 2019
accs_musfk01 2016, 2017, 2018, 2019
accs_mussm15 2013, 2014, 2015, 2016, 2017, 2018, 2019
accs_mustu17 2015, 2016, 2017, 2018, 2019
accs_mutsk02 2013, 2014, 2015, 2016, 2017, 2018, 2019
accs_mutsk09 2013, 2014, 2015, 2016, 2017, 2018, 2019
cik_Ben Courtney Creek 2017
cik_Big Creek 2014, 2015
cik_Copper River 2017, 2018
cik_Gibraltar River 2017, 2018
cik_Napotoli Creek 2014, 2015, 2016, 2017, 2018, 2019
cik_Neilson Creek 2014, 2015
cik_Panaruqak Creek 2015, 2017
cik_Roadhouse Creek 2018
cik_Silver Salmon Creek 2014, 2015, 2016, 2017, 2018
cik_Sivanguq Creek 2014, 2015, 2016, 2017
cik_Squaw Creek 2014, 2015
cik_Tunravik Creek 2014
cik_Yellow Creek 2015, 2016, 2017, 2018, 2019
fws_580223156504200 2015, 2016, 2017, 2018
fws_GELO 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2011, 2012, 2013
fws_KULR 2002, 2003, 2006, 2007, 2008, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018
fws_MALR 2003, 2006
fws_NCLO 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019
fws_Newhalen River 2009, 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018
fws_OSLR 2002, 2003, 2005
fws_OSLR2 2010, 2011, 2012, 2013, 2017, 2018
fws_PULO 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2012, 2014, 2015, 2016, 2017, 2018, 2019
fws_PULR 2003, 2004, 2005, 2008, 2009, 2010, 2011, 2012, 2013, 2014, 2015, 2019
fws_SLLR 2005, 2006, 2007, 2008, 2011, 2014, 2015, 2016, 2017, 2018
fws_TOLO 2003, 2004, 2005, 2006, 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018
fws_TOLR 2003, 2004, 2005, 2006, 2007, 2010, 2011, 2012, 2013, 2014, 2015, 2016
fws_WELR 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010
fws_WELR2 2012, 2013, 2014, 2015, 2016, 2017
NPS_Chulitna_SiteB 2008, 2010, 2011
nps_KATM_idavc_stream_water 2017, 2018, 2019
nps_KATM_lbrooo 2006, 2007, 2008, 2012, 2013, 2018
nps_KATM_margc_stream_water 2015, 2017, 2019
nps_KATM_naknlo 2006, 2007, 2009, 2010, 2011
nps_KATM_savor_stream_water 2018, 2019
nps_KATM_upatc_stream_water 2015, 2016, 2017, 2018, 2019
nps_LACL_kijilo 2016
nps_LACL_tazir_stream_water 2015, 2018, 2019
nps_LACL_tlikr_stream_water 2015, 2016, 2017
NPS_Little_Kijik_River_Dan 2001, 2002, 2003, 2004
NPS_Newhalen_River 2008, 2009, 2010, 2014, 2015, 2016
NPS_Six_Mile_Outlet 2007, 2008
usgs_15298040 2010, 2011, 2013, 2014, 2015, 2016, 2017
usgs_15300100 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013, 2014, 2015, 2016
usgs_15300250 2006, 2007, 2008, 2010, 2011, 2012, 2013, 2014, 2015, 2016
usgs_15300270 2012, 2013
usgs_15300300 2014, 2015, 2016, 2017, 2018, 2019
usgs_15300320 2011, 2012, 2013
usgs_15301500 2012, 2013, 2014
usgs_15302000 2014, 2015, 2016, 2017, 2018, 2019
usgs_15302200 2007, 2008, 2011, 2012, 2013, 2014, 2015, 2018, 2019
usgs_15302250 2007, 2008, 2009, 2010, 2011, 2012, 2013, 2014, 2015, 2018, 2019
usgs_15302300 2014
usgs_15302320 2014
usgs_15302812 2017, 2018, 2019
UW_Agulowak River 2017
UW_Agulukpak River 2011, 2012, 2016, 2019
UW_Aleknagik Bear Creek 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019
UW_Aleknagik Big Whitefish Creek 2011, 2012, 2013, 2014, 2015, 2016, 2018, 2019
UW_Aleknagik Eagle Creek 2011, 2012, 2014, 2015, 2016, 2017, 2018, 2019
UW_Aleknagik Hansen Creek 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018
UW_Aleknagik Happy Creek 2011, 2012, 2015, 2016, 2018, 2019
UW_Aleknagik Ice Creek 2011, 2012, 2017, 2018, 2019
UW_Aleknagik Mission Creek 2011, 2012, 2016, 2018, 2019
UW_Aleknagik Pfifer Creek 2011, 2012, 2016, 2017, 2018, 2019
UW_Aleknagik Silver Salmon Creek 2011, 2012, 2015, 2016, 2017, 2018, 2019
UW_Aleknagik Squaw Creek 2016
UW_Aleknagik Sunshine Creek 2014
UW_Aleknagik Yako Creek 2011, 2012, 2015, 2016, 2017, 2018, 2019
UW_Aleknagik Youth Creek 2016
UW_Beverley Hope Creek 2013, 2016, 2017, 2018
UW_Beverley Moose Creek 2013, 2015, 2016, 2017, 2018
UW_Beverley Silverhorn Creek 2013
UW_Beverley Uno Creek 2013, 2016
UW_Kulik Grant River 2016, 2017, 2018
UW_Kulik K-3 Creek 2012, 2013, 2015, 2018, 2019
UW_Kulik Kulik Creek 2016, 2019
UW_Little Togiak C Creek 2011, 2012, 2014, 2015, 2016, 2017, 2018, 2019
UW_Little Togiak Creek 2016, 2017, 2018, 2019
UW_Nerka Allah Creek 2011, 2012, 2014, 2015, 2016, 2017, 2018
UW_Nerka Bear Creek 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019
UW_Nerka Beaver Creek 2011, 2012
UW_Nerka Berm Creek 2012, 2013, 2014, 2015, 2016, 2018, 2019
UW_Nerka Bug Creek 2011
UW_Nerka Cabin Creek 2012, 2013, 2014, 2016, 2018
UW_Nerka Chamee Creek 2016
UW_Nerka Cottonwood Creek 2012, 2013, 2014, 2016, 2019
UW_Nerka Elva Creek 2013, 2014, 2016, 2017, 2018, 2019
UW_Nerka Fenno Creek 2011, 2012, 2016, 2017, 2018, 2019
UW_Nerka Hidden Lake Creek 2011, 2012, 2013, 2015, 2016, 2017, 2018, 2019
UW_Nerka Joe Creek 2011, 2012, 2013, 2014, 2015, 2016, 2018
UW_Nerka Kema Creek 2011, 2013, 2014, 2015, 2016, 2017, 2018, 2019
UW_Nerka Little Togiak River 2010, 2012
UW_Nerka Lynx Creek 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019
UW_Nerka Lynx Creek Cold Tributary 2011, 2012, 2013, 2014, 2016, 2017, 2019
UW_Nerka Lynx Lake Tributary 2017, 2019
UW_Nerka N4 Creek 2016, 2017, 2018, 2019
UW_Nerka Pick Creek 2012, 2013, 2014, 2015, 2016, 2018, 2019
UW_Nerka Rainbow Creek 2012, 2013, 2014, 2016, 2017, 2018
UW_Nerka Sam Creek 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2019
UW_Nerka Seventh Creek 2011, 2012, 2013, 2015, 2016, 2017
UW_Nerka Sixth Creek 2011, 2012, 2013
UW_Nerka Stovall Creek 2011, 2012, 2013, 2014, 2016, 2017, 2018, 2019
UW_Nerka Teal Creek 2012, 2013, 2014, 2015, 2016, 2018, 2019
mod_dat %>% 
  filter(Region == region_in) %>% 
  select(TempSens, log_slope, log_cat_elev, log_cat_slope, 
         wtd_slope_MEAN, snow_ind, wtd_lcld_jd, summer_precip) %>% 
  pairs(upper.panel = panel.cor, diag.panel = panel.hist, lower.panel = panel.smooth)

mod_dat %>% 
  filter(Region == region_in) %>% 
  select(TempSens, log_area, dist_coast_km, wtd_north_per,
         asrt_wet, asrt_glac, asrt_lake, glac_10) %>% 
  pairs(upper.panel = panel.cor, diag.panel = panel.hist, lower.panel = panel.smooth)

9 Prediction scenarios

Spatial domain: salmon habitat for all regions. Select HU12 that intersect anadromous waters catalog. Spatial scale: suggest using stream outlets associated with HU12 boundaries (smallest are 25 km2). We can use these outlets to create watersheds and calculate all the same covariates for prediction.

Temporal scenarios: Ideas for our scenarios include picking years within our temporal domain (2001-2019) that represent contrasting conditions that affect TS (e.g. low and high snow years). One problem with this idea is that conditions could vary across our study area so a high snow year in Bristol Bay may not match Copper. Alternatively, we could develop scenarios for prediction using high and low ranges of snowpack or precipitation for each region.

Comparison of standardized (subtract mean and divide by sd) temperature and precipitation values by different regions:

  • 2019 warm and low precipitation
  • 2016 warm and high precipitation
  • 2006 cold and high precipitation
  • 2017 and 2009 about normal years for both
clim_dat <- readRDS("output/clim_dat.rds")

#deviations from normal
clim_dat %>% 
  group_by(Region, parameter) %>% 
  mutate(Value_metric_sc = scale(Value_metric)) %>% 
  ggplot(aes(x = Year, y = Value_metric_sc, color = Region)) +
  geom_line() +
  scale_x_continuous(breaks = c(seq(2000, 2019, 2))) +
  facet_wrap(~parameter, scales = "free", ncol = 1) +
  geom_hline(aes(yintercept = 0)) +
  theme_bw() +
  theme(legend.position = "right") +
  labs(y = "Standardized values")

Snowtel data indicate April 1st SWE (magnitude of spring snowpack) for different sites across some of our regions. Note that there were no snowtel sites with data for Bristol Bay or Kodiak. From Tim’s paper, 2012 and 2013 were high snow years and 2015 was a low snow year.

  • 2012 was a high snow year
  • 2015 was a low snow year
snowtel <- readRDS("output/snowtel.rds")

#deviations from normal
snowtel %>% 
  filter(!grepl("Gulk|Tela", Snowtel_site)) %>% 
  group_by(Snowtel_site) %>% 
  mutate(swe_sc = scale(Apr1_SWE)) %>% 
  ggplot(aes(x = Year, y = swe_sc, color = Snowtel_site)) +
  geom_line() +
  scale_x_continuous(breaks = c(seq(2000, 2019, 2))) +
  theme_bw() +
  theme(legend.position = "right") +
  geom_hline(aes(yintercept = 0)) +
  facet_wrap(~Region, ncol = 1) +
  labs(y = "Standardized values")